arXiv:l504.05194v3 [astro-ph.SR] 10 Aug 2016 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 12 August 2016 (MN style file v2.2) 


Thermal Runaway During the Evolution of ONeMg Cores 
towards Accretion-Induced Collapse 

Josiah Schwab^’^, Eliot Quataert^’^, Lars Bildsten^’^ 

^Physics Department, University of California, Berkeley, CA 9J^720, USA 

^Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720, USA 
^Department of Physics, University of California, Santa Barbara, CA 93106 
"^Kavli Institute for Theoretical Physics, Santa Barbara, CA 93106 


Revision b07f982950elc5el44321cabb4d3e80735c26383 


ABSTRACT 

We study the evolution of degenerate electron cores primarily composed of the carbon 
burning products ^°Ne, and ^'^Mg (hereafter ONeMg cores) that are undergoing 
compression. Electron capture reactions on A = 20 and A = 2A isotopes reduce the 
electron fraction and heat the core. We develop and use a new capability of the Modules 
for Experiments in Stellar Astrophysics (MESA) stellar evolution code that provides 
a highly accurate implementation of these key reactions. These new accurate rates 
and the ability of MESA to perform extremely small spatial zoning demonstrates a 
thermal runaway in the core triggered by the temperature and density sensitivity 
of the ^°Ne electron capture reactions. Both analytics and numerics show that this 
thermal runaway does not trigger core convection, but rather leads to a centrally 
concentrated (r < km) thermal runaway that will subsequently launch an oxygen 
deflagration wave from the center of the star. We use MESA to perform a parameter 
study that quantifies the influence of the ^^Mg mass fraction, the central temperature, 
the compression rate, and uncertainties in the electron capture reaction rates on the 
ONeMg core evolution. This allows us to establish a lower limit on the central density 
at which the oxygen deflagration wave initiates of pc ^ 8.5 x 10®gcm“^. Based on 
previous work and order-of-magnitude calculations, we expect objects which ignite 
oxygen at or above these densities to collapse and form a neutron star. Calculations 
such as these are an important step in producing more realistic progenitor models for 
studies of the signature of accretion-induced collapse. 

Key words: white dwarfs - stars:evolution 


1 INTRODUCTION 


In this paper, we study the evolution of degenerate electron 
cores primarily composed of the carbon burning products 
isq, ^^Mg which are undergoing compression. 

Such objects can arise in several contexts: the late stages of 
evolution for super asymptotic giant branch (SAGB) stars 
(e.g.|Miyaji fc Nomoto||1987[ |Ritossa et al.|1999| iTakahasld] 


et al. 2013 Jones et al. 20131, where the compression is 


caused by the deposition of material from exterior shell¬ 
burning; in a binary system with a massive ONeMg white 
dwarf (WD) (e.g. Nomoto fc Kond^|1991 |, where the com¬ 
pression is caused by accretion from a non-degenerate com¬ 
panion; or as the remnant of a WD-WD merger, where the 
compression is caused by the cooling of the outer layers (e.g. 
Saio fc Nomoto|1985 ). 

As the core is compressed, the electron Fermi energy 
rises, eventually triggering exothermic electron capture reac¬ 


tions. Typically, exothermic captures on ^°Ne release enough 
energy to cause thermonuclear ignition of and formation 
of a deflagration. The final fate of the core (either explosion 
or collapse) is determined by a competition between the en¬ 
ergy release from the outgoing oxygen deflagration and the 
energy losses and decline in the electron fraction due to elec¬ 
tron captures on the post-deflagration material, which has 
burned to nuclear statistical equilibrium (NSE). The evolu¬ 
tion of these cores has been a subject of considerable previ- 


ous study (e.g. Miyaji et al.| 1980| Nomoto|1984| Isern et al. 

1991 

Canal et al.|1992||Gutierrez et al.|1996||Gutierrez et al. 

2005 

Jones et al.| 

20141). 


However, we revisit this topic (i) to test the effect of us¬ 
ing the state-of-the-art Modules for Experiments in Stellar 


Astrophysics (MESA) stellar evolution code ( Paxton et al. 
2011 2013| 20151, (ii) to demonstrate the effects of using the 


latest nuclear reaction rates |Martmez-Pinedo et al.|[2014|. 
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2 Schwab et al. 


(iii) to perform a more detailed parameter study of the ef¬ 
fects of a number of quantities, including the accretion rate 
M, magnesium mass fraction -^Mg, and initial core tem¬ 
perature, Tc, and (iv) to provide analytic estimates of the 
evolution up-to and including the onset of the oxygen defla¬ 
gration. 

In the present paper, we follow the common treat¬ 
ment in the literature and parameterize the evolution of 
ONeMg WDs as they approach the Chandrasekhar mass 
via compression of the outer layers. In future work we 
will assess whether the revised evolutionary model of WD 


merger remnants proposed by Shen et al. (20121 and Schwab 
|et al.| ( |2012[ ) modifies the likelihood of AIC in super- 
Chandrasekhar WD mergers. In § [^ we describe the treat¬ 
ment of weak reactions in the MESA code. In §|^we provide 
analytic estimates relevant to the evolution of the core. In 
§ l^we discuss the inputs to our MESA calculations and in 
§ present the results of these numerical simulations. § 
discusses the final fate of these cores. In § [^ we draw our 
conclusions and describe some important avenues for future 
work. 



2 WEAK REACTIONS IN MESA 


Weak reactions, specifically electron-capture and beta- 
decay, are central to the evolution of accreting degenerate 
ONeMg cores. The reduction in electron fraction (and cor¬ 
responding reduction in pressure) due to electron captures 
accelerates the contraction of the cores and the entropy gen¬ 
eration from these electron captures can directly ignite ther¬ 
monuclear reactions. 

This study makes use of MESA, a state-of-the-art open 


source code for stellar evolution calculations (Paxton et al. 


2011 20131. In particular we use the capability to calculate 


weak reaction rates directly from nuclear level and transition 
data, which is documented in the upcoming MESA Instru¬ 


ment Paper III (Paxton et al. 20151. This section summa¬ 


rizes the input data to this capability. The precise expres¬ 
sions which are evaluated as part of MESA’s on-the-fly weak 
reaction treatment are given in Appendix [A] 

We restrict ourselves to considering only a small set of 
A = 24 isotopes (^^Mg, ^^Na, ^'*Ne) and A = 20 isotopes 
(^°Ne, ^“O). Over the range of thermodynamic con¬ 

ditions encountered during the evolution of ONeMg cores, 
roughly 9 < logjQ p < 10 and 8 < logj^Q T < 9 (in cgs units), 
Takahara et al. (19891 identified the transitions that domi¬ 


nate the rate of each reaction. We consider only this limited 
set of transitions; they are listed in Table We have taken 
the comparative half-lives of these reactions from the up-to- 
date information compiled in Martmez-Pinedo et al. ( 2014[). 

In order to more easily visualize the data in Table |1[ 
we present energy level diagrams for the A = 24 (Fig. [r 
and A = 20 (Fi g. [2| nuclei. The se figures are modeled after 
those found in Takahara et al. (19891. The level structure 


of these nuclei is drawn from recent compilations of nuclear 


data (Tilley et al. 1998 Firestone 20071. We show all of 


the low-lying states that we consider, labeled with their J'^ 
(spin^'*'^'*^) values. The arrows indicate the limited set of 
transitions that we consider, which are only those which 
are “allowed” (Gamow-Teller: Ji = Jj,Jj ± 1, TTiiTj = 1; 
excluding Ji = J/ = 0). 


Figure 1. Energy level diagram for the A = 24 nuclei that we 
consider. The J'^ values are sometimes given an arbitrary offset 
(indicated via thin lines) in order to enhance legibility. The tran¬ 
sitions we consider are indicated with arrows. 
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Figure 2. Energy level diagram for the A = 20 nuclei that we 
consider. The values are sometimes given an arbitrary offset 
(indicated via thin lines) in order to enhance legibility. The tran¬ 
sitions we consider are indicated with arrows. 


3 ANALYTIC ESTIMATES 


Miyaji et al. (1980 \ provide a thorough discussion of the 


different phases of the evolution of an ONeMg core under¬ 
going compression. In order to gain some insight into the 
relevant physics, we first discuss a simple model of the evo¬ 
lution up until the onset of thermonuclear oxygen burning. 
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Evolution of ONeMg Cores 3 


Initial 

Final 

Qg 

Ei 

jr 

Ef 

TTT 

Jf 

loglo(/i) 


24 Na 

5.515 

0.000 

0.000 

0+ 

0+ 

0.472 

1.347 

1+ 

1+ 

4.815 

3.838 




0.000 

4+ 

3.972 

4+ 

6.209 


24 Ne 

2.467 

0.000 

4+ 

4.866 

3+ 

4.423 




0.472 

1+ 

0.000 

0+ 

4.829 




1.634 

2+ 

0.000 

2+ 

4.970 

20Ne 

20f 

7.025 

0.000 

0+ 

1.057 

1+ 

4.380 




0.000 

0+ 

0.000 

2+ 

9.801 

20p 

20 o 

3.815 

0.000 

2+ 

1.674 

2+ 

5.429 


1.057 

1+ 

0.000 

0+ 

4.211 


Table 1. The transitions used in the rate calculations. They are 
written as electron capture transitions, but the same transitions 
were used for beta-decay (swapping initial and final states). Qg is 
the energy difference between the ground states of the isotopes. 
E'l and E{ are respectively the excitation energies of the initial 
and final states, relative to the ground state, jr and J'^ are the 
spins and parities of the initial and final states. Allowed transi¬ 
tions do not have parity changes, (ft) is the comparative half-life 
in seconds, taken from [Martmez-Pinedo et al.| ( |M14[ l by dividing 
the constant 6144 s by their tabulated values of the transition ma¬ 
trix elements. The italicized (ft) value indicates an experimental 
upper limit; the effects of this transition will be discussed in § |5.3| 
All energies are in MeV. For level diagrams which illustrate the 
transitions, see Figs. and 


In discussing the analytic estimates below, we reference some 
of the numerical results from our fiducial MESA model for 
comparison. This model is a cold ONeMg WD (Ao = 0.5, 
Anb = 0.45, Aiug = 0.05) accreting at M = 10“® Mq yr“^. 


3.1 Overview of evolution 


We have a dense, degenerate core near the Chandrasekhar 
mass with a spatially-uniform composition of the carbon¬ 
burning products ^®0, ^°Ne and ^“^Mg, with mass fractions 
Ao, Anb, Aiug, respectively. Fiducially, we choose Ao = 0.5, 
Anb = 0.45, Xms = 0.05. This is similar to the central abun¬ 
dances observed in recent calculations of the evolution of 
intermediate mass stars that develop these cores (see e.g., 
figure 10 of Takahashi et al.||2013 |. Other recent models of 
super-AGB evolution ( Farmer et al.|20l5 1 show typical cen¬ 
tral magnesium fractions AMg ~ 0.03 in the cases where 
the carbon deflagration wave reaches the center (R. Farmer, 
private communication). 

The degenerate core is “accreting” at a rate M; such 
accretion might be set by carbon shell burning in an evolved 
star, accretion from a companion in a binary system, or 
cooling (and the concomitant reduction in pressure support) 
of the outer layers of a WD merger remnant. The key impact 
is that the core is being compressed on a timescale 


/din Pc 
\ dt 


f din pc 
VdlnM 


^ . ( 1 ) 
M 


For an object supported by degeneracy pressure and in hy¬ 
drostatic equilibrium, the central density rises rapidly as one 
approaches the Chandrasekhar mass. Therefore, the com¬ 
pression timescale is significantly shorter than the timescale 
for the growth of the core. For an ideal, zero-temperature 


white dwarf, in the range 9 < logj^Q pc < 10, 

~ 28 r V 

dluM \^10®gcm“3y 


( 2 ) 


which we obtained by calculating a sequence of models and 
fitting a power-law to the results. This implies 


^compress ~ 5xl0Vr 


10® g cm“® 


M 


10-®M© yr-i 


( 3 ) 


The dynamical time of the white dwarf is extremely short 

' 

tdyn ~ —f= « 10~^s 


CG-p 


10® g cm“® 


( 4 ) 


and so hydrostatic equilibrium will always be preserved (un¬ 
til collapse ensues, which we do not study in detail in this 
paper). 

The temperature of the core will be influenced by details 
of its previous evolution, such as the accretion history and by 
the abundances of isotopes which participate in Urea process 
cooling. However, if the compression timescale (and hence 
overall evolutionary timescale) is sufficiently slow, heating 
from compression and cooling from thermal neutrinos will 
reach a quasi-equilibrium ( |Paczyhski|1971| ) . Define the cool¬ 
ing time 


tcool - 


cpT 


( 5 ) 


where cp is the specific heat at constant pressure and Si, is 
the specific neutrino cooling rate. Then the relation tcooi = 
tcomprBss implicitly defines a temperature for a given density 
and will characterize the thermal state of the core aside from 
periods when e-captures rapidly release energy. In Fig.[^ we 
show this relation as a blue, dashed line and demonstrate 
that our MESA models (the black solid line) described in 
§[^and §1^ exhibit this relationship. 


3.2 Effects of electron captures 


As the core is compressed, the electron chemical potential 
increases. At zero temperature, the electron captures would 
occur when the Fermi energy reached the energy difference 
between the initial and final nuclear states. We refer to the 
density corresponding to this value of the chemical poten¬ 
tial as the threshold density; the terms sub-threshold and 
super-threshold reference this density. At non-zero temper¬ 
ature, even when the chemical potential is below this thresh¬ 
old, there are some electrons in the high energy tail of the 
Fermi-Dirac distribution which are available to capture. As 
a result, the electron capture rate has an exponential depen¬ 
dence on the density and temperature in the sub-threshold 
case. 


A simple form for the sub-threshold capture rate can 
be obtained by expanding equation (A15l in the limit that 
p + Q —kT (where p is the chemical potential and Q = 
Qg + Ei — Ef is the energy difference between the parent 
and daughter nuclear state) and assuming that the rate is 
dominated by a single transition that begins in the ground 
state. 


Abc 


21n2 / kT y 
(ft) VmeC®/ 



P + Q '\ 
kT ) 


( 6 ) 
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logio Pc [g cm 3] 


logio P [g cm 3] 


Figure 3. The black solid line shows the central density and tem¬ 
perature of the core as it is compressed with a surface accretion 
rate of M = 10 ~® MQyr~^ for approximately 20000 years of evo¬ 
lution. The red dashed lines indicate when the capture timescales 
for and ^^Ne become equal to the fiducial compression time 

of 10*^ yr. The blue dotted line shows where the neutrino cool¬ 
ing time and compression time are equal. This balance between 
compressional heating and neutrino cooling determines the ther¬ 
mal state of the contracting WD core (aside from brief periods 
when electron captures heat the core). The grey solid line shows 
where the energy generation from thermonuclear oxygen burning 
exceeds the thermal neutrino losses and we stop the calculation. 


Define the capture timescale to be the inverse of the 
electron capture rate ^capture = The onset of signifi¬ 

cant electron captures will occur when the capture time and 
the compression time become approximately equal. Setting 
^compress = ^capture gives an implicit relationship between p 
and T, which is a function of M. 

At zero temperature, the electron captures would occur 
at a density pec,o such that p(pec,o) + Q ~ 0. Solving equa¬ 
tion for pL and rewriting the solution in terms of p, we 
find that icompress ~ ^capture when 


Pec,0 


, 3fer 


In 2 In 2 


tc 


(ft) 


kT V f 


nieC^ 


kT 


(7) 


where we have neglected the much weaker density depen¬ 
dence of tcompress ItSelf. 

Equation 0 will be valid up until a temperature at 
which the transition rate from an excited state, suppressed 
by exp{—Ei/kT), becomes the dominant contribution to the 
rate. As a rule of thumb, for the transitions we consider, this 
will happen when T ~ Ei/{25k). 

Fig. i shows numerical solutions for the location in 
density-temperature space at which fcapture = 10 '* yr (which 
is approximately the compression timescale associated with 
an M = lO“®M 0 yr-*) for ^^Mg, ^^Na, ^^Ne, and ^“F. 


Figure 4. The solid lines show where the electron capture 
timescale is equal to 10 “^ yr (which is approximately the com¬ 
pression timescale of the WD core for M = lO“®M 0 yr“^). At 
p and T greater than those delineated by the solid lines, the 
capture time is less than the compression time. Each line is la¬ 
beled by the name of the isotope undergoing electron capture. 
The black dashed lines show the analytic approximation given in 
equation 0 - For ease of comparison with the analytic results, the 
Coulomb corrections discussed in Appendix are not present in 
these calculations. 


The approximations for the critical density based on equa¬ 
tion 0 are overlaid as dashed black lines and are in excel¬ 
lent agreement. The line for is always at lower density 
than that of ^'^Ne, meaning once the first capture in the 
^'^Ne —>■ ^°F —>■ chain occurs, the second will immedi¬ 
ately follow. This is not true for ^*Na relative to its parent 
^*Mg, meaning the captures in the ^*Mg —>■ ^*Na —>■ ^*Ne 
chain will happen at separate densities when log,Q T < 8.4. 

The electron captures also influence the temperature 
evolution of the core. When a capture occurs, the chemical 
potential of the captured electron, minus the change in nu¬ 
clear rest mass and the energy in the emitted neutrinos, is 
thermalized, heating the plasmaQ This heating is substan¬ 
tial, because the first capture is often into an excited state 
(meaning the chemical potential is higher when the rate of 
this transition becomes significant) and the second is typi¬ 
cally super-threshold. Does this heating drive convection? If 
so, this convection will efficiently transport the entropy out 
of the core while mixing in fresh fuel for electron captures. 

The electron captures generate entropy, creating a neg¬ 
ative radial entropy gradient in the core. The captures also 
reduce the electron fraction in the core, creating a positive 
radial gradient in V).. The entropy gradient is destabilizing, 
but the Ye gradient is stabilizing. Simulations which invoked 


the Schwarzschild criterion for convection (e.g. Miyaji et al. 


^ A more precise definition of the heating rate is given in Ap¬ 
pendix specifically equations | |A29[ | and ( |A30| |. 
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1980), which does not consider composition gradients, found 


that the captures do trigger convection. Simulations which 
invoke the Ledoux criterion, which does consider composi¬ 
tion gradients, found that the captures do not trigger con¬ 


vection (e.g. Miyaji & Nomoto 1987). Hence the different 


choices lead to qualitatively different evolution. 

The following calculation demonstrates why the elec¬ 
tron captures are unable to trigger convection when account¬ 
ing for stabilizing composition gradients. The Ledoux crite¬ 
rion for convective instability is 


Vad — Vt-|-B<0 


where 


B = -- 


Xt 


dlnP\ 

dlnYj 


P.T 


din He 
dlnP 


( 8 ) 


( 9 ) 


The captures occur over a narrow range in Fermi energy, and 
hence density. Therefore the gradients in T and Ve across the 
region where the captures occur will be large. This allows 
us to drop the Vad term. Replacing the gradients with finite 
differences, we then check the inequality 


A(lnT)>-— A (In He). 

XT VolnTe/pr 


( 10 ) 


For a cold plasma with degenerate electrons and ideal ions, 
{din P / dlnYe) ^ j, ~ 4/3 and xt ~ •ikT/{ZE-p)- If a mass 
fraction AX^ has undergone electron captures, the associ¬ 
ated change in temperature is 




( 11 ) 


where is the nuclear mass number of the species that 
is capturing and Ec is the average energy deposited by a 
capture. At the typical densities and temperatures in our 
calculation, the ions are a Coulomb liquid and so cp ~ 3k. 
The change in Y^ due to the captures is 

A7 

ATe RS — AA, . (12) 

-^C 

Both the A = 24 and A = 20 chains that we consider are 
two electron captures long, so we set AZ = —2. 


Substituting these estimates into equation (101 and sim¬ 


plifying, the condition for convective instability becomes 
E, 


El 


> 2 . 


(13) 


This inequality demonstrates that in order to trigger con¬ 
vective instability, the two captured electrons—which each 
have a characteristic energy of Ep —would have to deposit 
nearly all their energy as thermal energy. This is unrealistic, 
since substantial amounts of energy go into the rest mass of 
the daughter nucleus and to neutrinos. From our calculation 
of the heating rates, it is clear this inequality is far from be¬ 
ing violated: the A = 24 captures occur at Ep « 6.5 MeV 
and release E^ ~ 0.5 MeV; the A = 20 captures occur at 
Ep ~ 8.5 MeV and and release Ec ~ 3 MeV. Electron cap¬ 
tures do not directly trigger convection. 

A region which is Schwarzchild-unstable but Ledoux- 
stable is semiconvective. The semiconvective diffusion coef¬ 
ficient used in MESA ( [Paxton et al.|2013[ following [Langer] 
et al.|1983 1 is 

f K \ f Vt - Vad 


Dsc — Ois 


\6cpp 


)( 


B -|- Vad — Vt 


(14) 
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where K is the radiative conductivity. The values of Osc, the 
semiconvective efficiency adopted in the liter ature s pan the 
range 10“^ < Osc < 1 (Paxton et al. 

TM) [Yoon etaL][2^ . 


2013 


citing Langer 


Regions where the electron captures have not yet oc¬ 
curred and regions where they have completed do not have 
a Yc gradient. Therefore the width of the semiconvective 
zone Hsc will be roughly the length over which the density 
changes by an amount that shifts Ep by kT. We expect 


Ha, 




(15) 


where Hp is the pressure scale height. Defining / = Hgc/Hp, 
we find / ~ 0.03 in our MESA models, consistent with 
the above estimate. We define the timescale for semicon¬ 
vection to modify the composition and thermal structure in 
our models as 


tsc = -pT~ ~ 3 X 10 yr - 

-L^sc V Qisc 


m I ■ <“> 


For Oac ^ 1, isc in equation (16) is equal to or longer than 
time that elapses between ^^Mg captures and oxygen igni¬ 
tions in our fiducial model. Moreover, tsc is an upper limit: 
because of the thinness of the region with a YZ-gradient, an 
individual parcel spends less time in a semiconvective re¬ 
gion. Therefore, we do not consider semiconvection in our 
models. 

For realistic ^^Mg fractions, e-captures on ^"^Mg do not 
release enough energy to initiate thermonuclear fusion. As 
a result, the core continues to compress and we eventually 
reach the density where the captures begin on the A = 20 
nuclei. Once the capture on ^°Ne occurs, the capture on 
^°F occurs immediately. Like the A = 24 captures, the bulk 
of the energy deposition comes from this super-threshold 
capture, but in this case the energy per capture is sub¬ 
stantially greater, Ec ~ 3 MeV. The characteristic tempera¬ 
ture for oxygen ignition is approximately 10® K and so from 


equation (11), we estimate that oxygen will ignite after an 
amount AA^e ~ 0.1 has undergone capture. We halt our 
main MESA calculations when the energy generation rate 
from oxygen burning exceeds the cooling from neutrinos, 
implying that a nuclear runaway is assured. 

We have focused primarily on the evolution of the center 
of the core, but the density of the rest of the core increases 
during compression. Electron captures on the A = 24 ele¬ 
ments have been occurring off-center as parcels of the star 
reach conditions favorable for these captures. This is illus¬ 
trated in Fig. where one can see the depletion of ®^Mg in 
the inner O.3M0 of the star. 


3.3 Thermal runaway from ®®Ne Captures and 
Oxygen Deflagration Initiation 

In Fig. [^ we show the evolution of the center of our MESA 
models as the A = 20 captures beginj^The profiles are la¬ 
beled by the central heating time of the model, fheat.c = 
CpT/tnuc- At these temperatures, the energy generation rate 


^ The MESA run shown in this plot used a finer central spatial 
and temporal resolution than our fiducial case in order to better 
resolve the onset of these steep central gradients. 
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Figure 5. The structure of the model from Fig. [^at the end of the 
MESA calculation (when the energy generation rate from oxygen 
burning exceeded the neutrino cooling). The top panel shows the 
density (pg, p in units of 10^ g cm“^), the temperature (Tg, T in 
units of 10^ K), and the entropy per baryon s (in units of /c), as 
a function of enclosed mass. The bottom panel shows the mass 
fractions of ^'^Mg and ^^Ne as well as the electron fraction Ye. 
The small region in which there is a Ye gradient due to the A = 24 
captures has been moving outward in a Lagrangian sense. By the 
time the center reaches the density for ^®Ne captures, the inner 
O.3M0 has already been depleted of ^'^Mg due to e-captures.The 
subsequent evolution is discussed in §[^ 


Figure 6. The temperature (T) and electron fraction {Ye) pro¬ 
files in the MESA model shown in Fig.[^ as it is approaching the 
onset of O fusion. At these temperatures, the energy release is 
dominated by the A = 20 captures, so T and Ye are closely tied. 
The lines are labeled by the heating timescale at the center of the 
model. The dotted grey line in the top panel shows the slope of the 
temperature profile expected for a thermal runaway with diffu¬ 
sion. The mass resolution in this calculation is significantly higher 
than that in other figures in order to resolve the small region (reg¬ 
ulated by thermal diffusion) within which thermal runaway sets 
in. The dots on the orange (hottest) temperature curve indicate 
the locations of the innermost three MESA zones; the mass in the 
central zone is roughly 4 x 10“ 


is dominated by the A = 20 captures, which are undergo¬ 
ing a thermal runaway in the thermally conducting core. 


From equation (111, the change in Ye associated with in¬ 


creasing the temperature from its value before the A = 20 
captures, T 4 x 10® K, to the temperature for oxygen ig¬ 
nition, T ^ 10^ K, is AYe ~ 0.006, in good agreement with 
the change observed in the lower panel of Fig. Changes 
in T and Ye will no longer be so tightly coupled once energy 
release from oxygen fusion exceeds that from electron cap¬ 
tures, pushing the core towards convective instability. How¬ 
ever, in order to reach convective instability, equation (101 
requires 


AT > 

Sk 


3.5 X lO^'K 


AKe 

0.006 


(17) 


a temperature so large that the central heating timescale 
from oxygen fusion would be fheat.c ~ 10~® s. We show here 
that a thermal runway is triggered long before such a con¬ 
dition is reached. 

The A = 20 electron captures occur in an environment 
where the electron Fermi energy is below the energy thresh¬ 
old. In this sub-threshold case, those electrons that capture 
are on the thermal tail of the distribution, makii^ the rate 
very sensitive to both density and temperaturelj We now 


® Because the bulk of the heating comes from the super-threshold 


show that this naturally leads to a local thermal runaway in 
the core whose size is limited by thermal conduction. This 
runaway provides the “hot-spot” needed to initiate the oxy¬ 
gen deflagration from the center of the star. 

The strong density sensitivity of the A = 20 captures 
implies that the runaway will begin at the exact center of 
the isothermal core. However, the pressure declines away 
from the core, leading to a temperature gradient on the scale 
over which the electron capture rate (and hence the heat¬ 
ing rate) varies by order unity. In this sub-threshold case, 
din X/dln P — Ep/h^kT), so the change in pressure needed 
to have the rate be less at the outer edge than the center is 


AP Rs 4Pc 


kT \ ^ f PcYe 

EfJ~ 


kT, 


(18) 


Hydrostatic equilibrium provides such a pressure change 
over a length scale 


\2ttGpI 


(19) 


where we have used Tc « 4 x 10® K and p, ~ 9x 10® g cm 


electron capture on that immediately follows the ^®Ne cap¬ 
ture, the capture rate on ^'^Ne is a good proxy for the temperature 
and density dependence of the heating rate. 
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Evolution of ONeMg Cores 7 


corresponding to the onset of j 4 = 20 captures in our fidu¬ 
cial model. This estimate is consistent with the length scale 
observed at the onset of the runaway in Fig. 

The subsequent evolution of the runaway is driven by 
the temperature sensitivity of the sub-threshold electron 
capture rate. This rate is well-approximated by equation Q, 
yielding a logarithmic derivative of the rate with respect to 
temperature of 


dlnA_ v + Q 

d\nT ^ 


( 20 ) 


Physically, the second term is the how far the transition is 
below its threshold energy (in units of kT). As in the A = 
24 case, electron captures become important when ~ 
tcompress; from equatiou this occurs at (/i + Q)/{kT) « 
— 14. The thermal runaway is sufficiently rapid that v + Q 
remains approximately fixed. This implies that the captures 
will be extremely temperature sensitive, scaling as A oc T", 
where 


_ din A 
^ dlnT 


3-f 14 


T 


4 X 10« K 


( 21 ) 


In the following estimates, we will take n « 12, and since 
n S> 1, we will treat n « n ± 1. 

Hence, as captures begin, their density dependence 
leads to a temperature gradient on the length scale given 


by equation (19l. Because convection is not initiated, the 


heating remains local, and this temperature gradient will 
grow with time in a thermal runaway. Once it is sufficient 
to cause an order unity variation of the capture rate across 
a given length r, the gradient will become non-linear. The 
hotter part will begin to evolve more rapidly and the evo¬ 
lution of the cooler part will freeze-out. This will occur 
when dT/dr « {Tln)/r and so on its own, thermal runaway 
leads to a characteristic profile where d In T/d In r « 1/n. 
However, thermal conduction limits the volume that can 
runaway to a fixed temperature, keeping regions where 
Uh f^heat approximately isothermal. The thermal diffusiv- 
ity from electron conduction is Dth ~ 30 cm^ s“^(T/10® K), 
meaning that the timescale for conduction to modify the 
thermal structure over a lengthscale r is 

Therefore, the size of the isothermal region at the center 
of the model scales like r oc Thus, as the run¬ 

away progresses, it will create a temperature profile with 
din T/d In r ~ —1/5. The dotted grey line in the top panel 
of Fig.[^shows this slope, which agrees well with the temper¬ 
ature evolution in the MESA calculations. The semiconvec- 
tive instability grows on the thermal diffusion time. During 
the thermal runaway, by definition, fheat Uh- Hence, the 
evolution of the core during this phase will be sufficiently 
fast that semiconvection will not modify the temperature or 
composition. 

This thermal runaway leads to a small volume at the 
core reaching very high tempertures, eventually to values 
large enough for heating from oxygen fusion to play a role. 
This occurs when T « 1.1 x 10® K, where the heating time is 
fheat ~ 10“® yr. From equation (221, the hottest (isothermal, 


so tth ~ Iheat) part of the core will have a size r « 3 x 
10® cm, which encompasses about 3x 10“®® of the total mass. 


The finest central zoning that we were able to achieve in 
our MESA calculations (as shown in Fig. was a mass 
resolution of approximately 4 x 1O“®®M0. Therefore, just 
as oxygen burning begins to dominate the energy release, 
the small size of this region prevents us from continuing to 
follow its evolution in our full star MESA simulations. 

The conditions created in the core of the star as energy 
generation by oxygen fusion begins to dominate over A = 20 
captures lead naturally to the development of an oxygen de¬ 
flagration wave. In particular, we have shown that oxygen 
fusion begins in a region at the core of the star whose size 
is determined by tth < theat- With theat identified as the 
heating time associated with oxygen fusion, this is precisely 
the condition for the onset of an oxygen deflagration wave; 


Timmes & Woosley (19921 defined the deflagration “trigger 


mass” to be the mass contained within the region satisfy¬ 
ing this constraint. Therefore, we are confident that the hot 
central region present at the end of our MESA calculations, 
being unstable to thermal runaway, will continue to grow in 
temperature and shrink in siz e, eventually reaching t he 1am - 
inar deflagration solutions of Timmes & Woosley (1992 \ ^ 
The outgoing deflagration wave will sweep across this ther¬ 
mally unstable core in less than one second. 

It is important to stress that the onset of the oxygen de¬ 
flagration in the AIC context is substantially different than 
the “simmering phase” in single degenerate Type la super¬ 
novae progenitors. There, after pycnonuclear carbon ignition 
occurs, the entropy release from carbon burning drives the 
formation of a central convection zone. The growth and heat¬ 
ing of this convective zone lead to a significant decrease in 
the central density between the time of carbon ignition and 
the development of a deflagration. 

In our models, by contrast, no central mixing occurs 
because of the stabilizing effect of the composition gradient 
associated with A = 20 captures. Therefore the central den¬ 
sity at which oxygen ignition occurs, and at which we halt 
our MESA calculations, is a good estimate of the central 
density at which the oxygen deflagration develops. We dis¬ 
cuss the propagation of this deflagration and its influence 
on the final outcome in § [^ 


4 DETAILS OF MESA CALCULATIONS 

All of the calculations performed in this paper are based on 
revision 6596 (released 2014-06-07), with some modifications 
to support our weak rate calculations. The incorporation of 
these changes into the mainline MESA code will be discussed 
in the upcoming MESA Instrument Paper HI (Paxton et 
al. 2015). As required by the MESA manifesto, the inlists 
and source code modifications necessary to reproduce our 
calculations will be posted on http://mesastar.org 


4.1 Generation of Initial Models 

In order to perform the parameter study discussed § [^ it 
is necessary to have a set of models of ONeMg cores with 


^ At the density in our MESA calculations, the laminar deflagra¬ 
tion width is (5 Ri 3 X 10“® cm, far below our ability to resolve in 
our full star simulations. 
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a range of different temperatures and compositions to use 
as initial conditions. We generate an idealized set of models 
via the following ad hoc steps. During each step, all nuclear 
reactions are turned off, ensuring that the model will con¬ 
tinually contract until halted by degeneracy pressure. 

We begin with a 1.325Mq pre-main sequence model of 
normal (roughly solar) composition. We evolve this model 
until it reaches a central density of logjg P = 3 (cgs). We 
then relax the (homogeneous) composition to the desired 
isq, 20 ]\jg^ ^^Mg mass fractions and allow the model to 
evolve until the central density reaches logj^Q p = 7. Then we 
set the model to accrete at a constant M and evolve until 
the central density reaches logj^g p = 9.4, which is still below 
the threshold for the onset of the electron capture reactions 
of interest. In order to achieve different core temperatures, 
we vary M ; models with higher (lower) accretion rates have 
less (more) time for neutrino cooling to carry away energy 
and are correspondingly hotter (colder). By this means, we 
arrive a set of models with varied compositions and central 
temperatures to use as initial models. 


4.2 Important MESA Options 


While our full inlists will be made publicly available, we 
highlight some of the most important MESA options used in 
the calculations. This section assumes the reader is familiar 
with sp ecific MESA options. Plea se consult the instrument 

20131 and the MESA websit^ 


papers (Paxton et al. 


2011 


for a full explanation of the meaning of these options. 

Since MESA is an implicit code, it is important that we 
choose timesteps that will resolve the processes of interest. 
The evolution of the ONeMg cores is driven by the increase 
in central density (and hence Fermi energy) caused by the 
ongoing compression. Therefore, our default runs include a 
timestep criterion based specifically on changes in central 
density 


delta_lgRho_cntr_hard_limit = 3e-3 
delta_lgRho_cntr_limit = le-3 

in addition to the primary spatial and temporal convergence 
settings of 


varcontrol_target = le-3 
mesh_delta_coeff =1.0 . 


Evidence demonstrating that this set of MESA options 
yields a converged result is shown in Appendix [C] 

These calculations use a nuclear network based on the 
co_burn.net network included with MESA with the addi¬ 
tion of the isotopes ^^Ne, and ^'^Na and the weak 

reactions linking the A = 20 isotopes to ^°Ne and the A = 24 
isotopes to ^^Mg. The special treatment of these weak re¬ 
actions (as discussed in Appendix]^ is activated with the 
options 


use_special_weak_rates = .true. 
ion_coulomb_corrections = ’PCR2009’ 
electron_coulomb_corrections = ’Itoh2002’ 


where the last two lines select the Coulomb corrections dis¬ 
cussed in Appendix [B| 


® http://mesa.sourceforge.net 


5 PARAMETER STUDIES 


In this section, we use a suite of MESA calculations to study 
how a variety of parameters affect the evolution of these 
cores. The key question we will answer is whether reasonable 
variation in these paramet ers w ill affect the final outcome. 

The first parameter (§ 5.11, the initial ^'’Mg mass frac¬ 
tion (XMg), is an intrinsic property of the ONeMg core, set 
during the process that produced the core. Variation in this 
value may reflect variation in the formation process (e.g., the 
initial mass of the star that produced it) as well as limits of 
our knowledge (e.g., uncertainties in quantities such a s th e 
^^C(q, 7)^®0 reaction rate). The second parameter (§ 5.21, 


the accretion rate M, is set by the current state of the sys¬ 
tem (e.g., the properties of a binary companion, the details 


of shell-burning). The third parameter (§ 5.31, the strength 
of the second forbidden transition between the ground states 
of ^®Ne and ^°F, reflects a limit in our current knowledge. 


5.1 Effect of a ^"^Mg mass fraction 


Gutierrez et al. (20051 performed a parameter study of 


the effects of the ^"‘Mg mass fraction. We follow their ap¬ 
proach of varying the central ^^Mg fraction, while holding 
the ^®0/^°Ne ratio fixed[^We explore a wide range of ^""Mg 
mass fractions, from 0.01 up to 0.20. This latter value is 
well above the expected ^‘*Mg fraction given current reac¬ 
tion rates. Our results are shown in Fig. 

The temperature increase due to the A = 24 electron 
captures scales roughly linearly with XMg, as expected from 
equation (111. For XMg 0.07, neutrino cooling erases the 


effect of the heating and the trajectories converge back to¬ 
wards the tcooi = fcompress relation described in § [^ Cor¬ 
respondingly, the density at which the captures on ^°Ne 
occur—and thus the density at which oxygen ignites—is in¬ 
dependent of XMg. For XMg ^ 0.07, an increase in XMg 
leads to the onset of ^°Ne captures (and oxygen ignition) 
at a higher density. At even higher values (not shown), the 
heating from the A = 24 captures is sufficient to directly 


ignite oxygen burning as noted by Miyaji & Nomoto (19871 


and Gutierrez et al. (20051. 


In the limited set of models presented in [Gutierr^ 


et al. (20051 this bifurcation in the core evolution around 


Aivig ~ 0.07 is not evident. However, it has a clear physi¬ 
cal explanation. In Appendix [D| we discuss a simple model 
of a zero-temperature white dwarf with a low-Ve core that 
explains these results; here, we demonstrate the consistency 
of these calculations with our MESA models, eliding the de¬ 
tails. 

Fig. [8] shows the compression time in two models, one 
on each side of this threshold value of Aivig. The core evo¬ 
lution for both models follows the dashed black line de¬ 
fined by equation § Up until the onset of the A = 24 
electron captures at log^Q pc ~ 9.6. Above this density, the 
= 0.05 model (blue line) follows the dotted line, which 
is the track expected for the fiducial value of M and the 


° In accordance with our fiducial model, we set this ratio be at 
10/9. The choice of this ratio does slightly influence the density 
at which captures occur through the dependence of the Coulomb 
corrections on Z. However, as a small shift on top of a small shift, 
we do not explore variations in this ratio. 
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Figure 7. The evolution of the central density and temperature 
for different magnesium mass fractions. At lower JAjvig values, the 
density at which oxygen ignition occurs (the end of the track) is 
independent of A]y[g; at higher XMg values, the density at which 
oxygen ignition occurs increases with increasing AMg- The text 
and Figs. 1^ and [^explain the origin of this trend. 


value of d In M/d In Pc in equation Q calculated from a zero- 
temperature white dwarf model in which Ye decreases for 
logj^g p > 9.6. See Fig. |D1| and surrounding discussion for 
the details of this zero-temperature model. 

In the AMg = 0.15 model (yellow line in Fig. [^, once 
the ^'^Mg captures occur at logjQ Pc ~ 9.6, the compres¬ 
sion timescale begins to fall dramatically. By the time ^°Ne 
capture densities are reached (logjg p « 10), the compres¬ 
sion timescale is orders of magnitude smaller than in the 
lower Aiug models, though it remains significantly longer 
than the dynamical time. Recall that significant electron 
captures only occur when the capture time satisfies the re¬ 
lation tcompress = tcapture; this means that for shorter com¬ 
pression timescales, the core must reach higher densities, and 
hence higher capture rates, before the effects of the captures 
become apparent. As Aiug increases, models experience a 
larger drop in Ye, and compress more quickly. This explains 
the trend of increasing oxygen ignition density with increas¬ 
ing Aiug seen in Fig. 

To physically understand the different evolution of the 
Aiug ^ 0.07 models, we consider an idealized model of the 
effect of Mg captures on the structure and stability of the 
ONeMg core. We assume that the A = 24 electron captures 
occur instantaneously above a density of p„ = 9.6. Using the 
approach described in Appendix]^ we can then determine 
the central density of the zero temperature model with the 
maximum mass. The result of this calculation is shown as 
the dashed line in Fig. Any model with Aiug // 0.07 will 
cross the stability line before the onset of ^°Ne captures. 
Moreover, for a larger change in Ye (associated with a larger 
Aiug in the current example), the onset of instability occurs 


Figure 8. The compression time of selected MESA models. The 
model with XMg = 0.05 is compressed at a rate controlled by the 
accretion rate. The dashed line shows the compression rate given 
in equation § and the dotted line shows the compression rate 
expected for a zero temperature white dwarf in which Ye sud¬ 
denly falls at logj^Q p = 9.6, due to the electron captures on ^^Mg 
and ^^Na. The agreement demonstrates that the central density 
is controlled by the total mass. The model with AMg = 0.15 ex¬ 
periences a much more dramatic decrease in the compression time 
because of the larger decrease in Ye (see Fig. [^. The black dash 
dotted line shows the neutronization timescale expected from the 
calculations in Appendix]^ where the central density is evolving 
at fixed mass. Note that in both cases the compression timescale 
still remains orders of magnitude longer than the dynamical time. 


at lower central density. These results explain the qnalita- 
tively different behavior of the high Aivig models in Fig. 

Above the dashed line in Fig. the zero temperature 
models are dynamically unstable and would contract on the 
dynamical timescale. But the characteristic electron cap¬ 
ture timescales are longer than the dynamical time, and so 
the assumption that the captures are effectively instanta¬ 
neous (used in the idealized models in Appendix D) does not 
hold in the real MESA models. As the contraction timescale 
gets shorter, only material at densities where the capture 
timescale is shorter than the contraction timescale can have 
had signihcant captures. Therefore the density above which 
the captures have completed, pn, shifts to higher values. 
There is no longer time for the total mass to change and so 
the timescale for the evolution of the central density is no 
longer set by the accretion rate. Instead, the core compresses 
on the significantly shorter neutronization timescale. 


in — 


fdlnYe 
\ dt 


FeA; 


Mg 


Aiug A, 


80 


Mg A ec 


Ai 


Mg 


0.15 


t. 


capture • 


(23) 


Given a fixed M, the models in Appendix give a rela¬ 
tionship between pe and Calculating in by evaluating 
fcapture at the pn Corresponding to each pc gives the dash 
dotted line in Fig.[^ which agrees well with the result of the 
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Figure 9. The evolution of the central electron fraction and cen¬ 
tral density for two values of X]vig (blue and yellow solid lines). 
The time elapsed during the evolution from log^Q p = 9.65 to the 
onset of ^®Ne captures is indicated next to each track. The dashed 
line shows the stability curve for a zero temperature white dwarf 
which neutronizes to the value of Ye,c shown on the x-axis at a 
density of log^^Q p = 9.6. The dotted red lines show the threshold 
electron capture densities for ^®Ne and ^'^Mg. For Xjvig ^ 0.07, 
the onset of A = 24 electron captures reduces Ye such that subse¬ 
quent compression drives the equivalent zero temperature mod¬ 
els into a dynamically unstable region of parameter space. Past 
this point, the contraction accelerates significantly (as shown in 
Fig.||. In these cases, captures alone have assured collapse. 


MESA calculation. See Fig. |D2| and surrounding discussion 
for the details of these zero-temperature models. 


5.2 Effect of central temperature and accretion 
rate 

As illustrated in Fig. the density at which electron cap¬ 
tures begin is temperature dependent. Our fiducial model 
begins at a central density logjQ pc ~ 9.4 and log,^Q Tc ~ 8.4. 
This central temperature is a free parameter, but as dis¬ 
cussed in § [^ a new central temperature will be established 
by the balance between neutrino cooling and compressional 
heating. Therefore, the central temperature when captures 
occur (in particular the A = 20 captures, and quickly there¬ 
after oxygen ignition) is weakly dependent on the initial tem¬ 
perature. This fact makes it difficult to separately illustrate 
the effects of the temperature and accretion rate. 

In Fig. we show the evolution of our fiducial model 
with 4 different accretion rates. The onset of captures is 
less temperature sensitive than one would infer from Fig. 
This is because at a higher M, while the quasi-equilibrium 
core temperature is higher (increasing the electron capture 
rates), the compression time is also shorter, and so the den¬ 
sity at which fcompress ~ fcapture Buds Up having a weaker 
dependence on the accretion rate. At the lowest accretion 
rate shown in Fig. 10 (M = lO~®M0yr“^), the evolution 


Figure 10. The fiducial model evolved with different values of 
M. The central density at which the A = 20 captures occur de¬ 
pends weakly on the accretion rate. At M = lO“®M 0 yr“^ the 
central temperature remains low enough that the ^^Mg and ^^Na 
captures occur at two separate critical densities. The dependence 
of the oxygen ignition density on M is weak. 


appears qualitatively different. Looking at Fig.[^ this is be¬ 
cause the central temperature remains sufficiently low that 
electron captures on ^"^Na do not occur immediately after 
electron captures on ^'‘Mg, but are delayed until higher den¬ 
sities (logio Pc ~ 9Y). 

Fig. demonstrates the independence of the oxygen 
ignition density on the initial central temperature. These 
models begin right before the A = 24 captures, at logiQ pc ~ 
9.55, so that the core temperature does not change sub¬ 
stantially before the onset of the captures. By the time the 
A = 20 captures occur, the temperature differences have 
been erased by neutrino cooling and compressional heating, 
as the core evolves towards the tcompress = fcooi thermal 
state discussed in §[^ As a result, there is little effect on the 
density at which oxygen ignition occurs. 


5.3 Effect of a ^^Ne forbidden transition 


Martmez-Pinedo et al. (20141 discuss the non-unique second 
forbidden transition from the 0^ ground state of ^°Ne to 
the 2"*" ground state of ^°F. The matrix element for this 
transition only has an experimental upper limit. They show 
that this transition can potentially dominate the rate for 
temperatures less than 9 x 10* KIj 


^ The results of both [Martmez-Pinedo et al. | |2014| | and of this 
work are obtained by treating the phase space factor of this second 
forbidden transition as that of an allowed transition. As discussed 
by |Martmez-Pinedo et ak] | |2014| |, the true shape factor could con¬ 
tain additional powers of the energy, which would further increase 
the rate, and can potentially offset the possibility that the matrix 
element is below the current experimental upper limit. 
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Figure 11. An illustration that the oxygen ignition density is 
independent of the initial central temperature. The models begin 
with log^o pc ~ 9.55, before the A = 24 e-captures. The grey 
dashed lines show when ^capture = ^compress for ^'^Mg and the 
models show the expected temperature dependence for the A = 24 
captures. However, by the time the A = 20 captures occur, the 
temperature differences have been erased by neutrino cooling and 
compressional heating, and thus there is little effect on the density 
at which oxygen ignition occurs. 


CM 



logio Pc [g cm 3] 


Figure 12. The effect of the second forbidden transition from 
the 0+ ground state of ^'^Ne to the 2+ ground state of The 
dotted, dashed, and dash-dotted lines show where the timescale 
for ^'^Ne captures is equal to the fiducial compression timescale 
(10^ yr) for different values of the matrix element (see plot leg¬ 
end) . The solid lines of matching color show the evolution of the 
fiducial MESA model using these rates. While the onset of ^®Ne 
captures shifts significantly if the O'*" —> 2+ transition is at the 
experimental upper limit, the shift in the density at which oxygen 
ignition occurs is substantially smaller. 


This transition can affect the critical de nsity at which 
'^^’Ne captures begin. The broken lines in Fig. |l2| shows the 
critical curves for '^''Ne capture obtained by setting the cap¬ 
ture rate equal to the fiducial compression rate, correspond¬ 
ing to M = 10“® Mq yr”*^. With the matrix element at the 
current experimental upper limit, the onset of captures is 
shifted to lower density (0.15 dex in logj^pp). At a value a 
factor of below the upper limit, the shift is very approx¬ 
imately halved (depending on the temperature). At a value 
a factor of 10® below the upper limit, the transition ceases 
to have a substantial effect. 

The solid lines in Fig. show the evolution of our 
fiducial model with each of these different choices for the 
strength of this transition. While the choice substantially af¬ 
fects the onset of ■^®Ne captures, it has a less significant effect 
on the density for oxygen ignition. Unlike the other transi¬ 
tions, which reach the critical capture timescale while they 
are sub-threshold, this transition is super-threshold. Cor¬ 
respondingly, the electron capture rate is less temperature 
sensitive. Its less rapid increase, coupled with the compres¬ 
sion timescale dropping due to the decrease in Ye, gives time 
for the core density to increase before the onset of oxygen 
ignition. 

In the calculation with the transition at the upper limit 
(solid yellow line), the central temperature does not reach 
the oxygen ignition line. This is because the thermal struc¬ 
ture of the remnant is such that ignition occurs mildly off- 
center. Our future calculations will determine whether this 
has any effect on the ensuing evolution. 


6 DISCUSSION 

As described by |Miyaji et al.| | |1980[ ), the final outcome of 
an ONeMg core as it approaches the Chandrasekhar mass, 
either explosion or collapse, is determined by a competition 
between the energy release from the outgoing oxygen defla¬ 
gration and the energy losses due to electron captures on 
the post-deflagration material, which has burned to nuclear 
statistical equilibrium (NSE). 

As discussed in § |^ the small length scale of the defla¬ 
gration means that we are unable to follow this phase with 
the MESA calculations presented in this paper. However, in 
lieu of a full calculation, we present a few order-of-magnitude 
estimates relevant to the outcome. 

At the time of collapse, the total energy of our fidu¬ 
cial white dwarf is E ~ —6 x 10®® erg. Oxygen burning to 
NSE yields approximately 1 MeV per baryon, meaning the 
energy release from burning 0.3 Mq of material can unbind 
the white dwarf. This energy release is required for the defla¬ 
gration wave to significantly change the structure of the star. 
Prior to the deflagration wave burning through ~ O.SMq of 
material, the structure of the WD core will remain relatively 
unchanged unless electron captures cause collapse. If the de¬ 
flagration moves at some fraction / of the sound speed, the 
timescale for it to propagate though the central O.SMq is 


td 



dr 

fCs 



(24) 
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where we have evaluated the integral using the structure of 
our MESA model at the end of the calculation. 


Nomoto & Kondo (19911 found that the critical de¬ 


flagration speed that demarcated the boundary between 
a model that explodes and a model that collapses was 
/ « 0.03. The work of Timmes & Woosley (19921, which 
simulated conductively-propagating oxygen deflagrations in 
detail gives a htting formula for the laminar deflagration 
speed of an oxygen flame of 


Vd = 51.8kms 


6 X 10® g cm~® 


0.6 ) 


(25) 


At p « 9 X 10®gcm“® and Xq = 0.5, this gives Vd ~ 
70kms“^, which corresponds to / = Vd/cs ~ 0.005. Based 
on an analysis of the growth of the Rayleigh-Taylor insta¬ 
bility, they conclude that these conductive flames are likely 
to remain stable. Therefore the laminar flame velocity is 
representative of the true flame speed in the inner part of 
the star. In particular, see figure 10 in [Timmes fc Woosley | 
(19921, noting that R{Mr = 0.3Mq) ~ 300 km. 

Based on these flame calculations, as well as several KE¬ 
PLER simulations using these speeds, [Timmes fc Woosleyl 


(19921 concluded that above a core density of 9 x 10® gem ^ 


the white dwarfs should collapse to a neutron star. The 
lowest central density at which oxygen ignition occurred 
in our parameter study (§ was logj^Q pc = 9.93. That is 
pc ~ 8.5 X 10®gcm“^, which is only marginally below this 
critical value. 

The timescale on which the core is neutronizing due to 
electron captures on the NSE-composition material can be 
written as tn = {dhiYc/dt)~^. The methods presented in 
this paper are not appropriate for calculating weak rates in 
NSE material. Instead, we take Ye from tables generated 


by Seitenzahl et al. (20091. Fig. 13 shows the neutroniza- 


tion timescale as a function of density and temperature for 
Ye = 0.49, the approximate central value in our fiducial 
model at oxygen ignition. Once the deflagration forms, the 
density of the post-deflagration material is less than the cold, 
upstream material. The MESA models reach oxygen igni¬ 
tion at logiQ Pc ~ 10, where this density change is small, 
Apjp « 0.1 ( Timmes fc Woosley[|1992 1. Therefore, the den¬ 
sity of the post-deflagration ash will be approximately the 
same as the density at which oxygen ignites, so Fig. [^indi¬ 
cates that the relevant neutronization timescale is approxi¬ 
mately 0.2 s. 

This estimate of the neutronization timescale is suffi¬ 
ciently shorter than the timescale on which the deflagration 
wave unbinds the star (equation [24[ | that it suggests that the 
end result of oxygen ignition following e-captures on ®®Ne 
will be collapse to a NS rather than a thermonuclear explo¬ 
sion. Future work will clarify this in the context of full-star 
simulations. 


7 CONCLUSIONS 


We have provided an updated analytic and numerical un¬ 
derstanding of the evolution of accreting and compressing 
ONeMg cores up to the initiation of oxygen burning in the 
core. This study was enabled by new capabilities of the 
MESA (Paxton et al. [201 1| 2013 20151 stellar evolution 
code. In particular, we have implemented a highly accurate 



Figure 13. The neutronization timescale for Ye = 0.49, which 
is roughly the central value of Ye at the end of our fiducial cal¬ 
culation. Contours are labeled by timescale. This uses the NSE 
electron capture rates from [Seitenzahl et ah] ( [2009^ . Since our 
compressing ONeMg cores reach oxygen ignition at logj^Q pc ~ 10, 
the relevant neutronization timescale is approximately 0.2 s. This 
is less than the timescale for the O deflagration wave to release 
enough energy to unbind the star (equation [24[ l, suggesting that 
collapse to a NS will ensue. 


treatment of the key electron capture rates (on A = 20 and 
24 nuclei) using modern microphysics from|Martfnez-Pinedo 


et al. (20141, which is summarized in Appendix |A| 


We have demonstrated analytically and numerically 
that neither ^''Mg or ®®Ne captures release sufficient heat 
to generate convection in the core. Instead, the core under¬ 
goes a thermal runaway triggered by the energy released by 
®®Ne captures. This centrally concentrated runaway initi¬ 
ates oxygen burning and launches an outgoing oxygen defla¬ 
gration wave at a time when the central density is at least 
8.5 X 10® g cm“®. Based on order of magnitude estimates and 
previous work of Timmes & Woosley (19921, we expect ob¬ 
jects which ignite oxygen at such high densities will collapse 
and form a neutron star due to continued electron captures 
on the NSE ashes produced by oxygen burning. 

Given the sensitivity of the final outcome of compress¬ 
ing ONeMg cores to the central density at the time oxy¬ 
gen burning begins (see we also performed a parame¬ 
ter study which demonstrated the influence of a number of 
factors on this density. We inv estig ated the effects of vary¬ 
ing the initial ®"^Mg fraction (§ 5.11, the initial central tem¬ 
perature and accretion rate (§ 5.21, as well as the poten¬ 
tial inclusion of a particular forbidden transition (§ 5.31. 
Figures [^ and demonstrate that values of AMg 0.07 
cause the core to contract more rapidly after A = 24 the 
captures, which leads to oxygen ignition at higher densities 
(thus further favoring collapse to form a neutron star). We 
also demonstrated the importance of the balance between 
neutrino cooling and compressional heating in setting the 
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central temperature of ONeMg cores during much of their 
evolution | Paczyhski|1971 1. This implies that the core typ¬ 
ically loses memory of its initial central temperature. 

The approach of an ONeMg core to the Chandrasekhar 
mass is relevant to the late stages of evolution for super 
asymptotic giant branch stars, binary systems with an ac¬ 
creting ONeMg WD, and the remnants of WD-WD mergers. 
Our calculations here are an important step in producing 
more realistic progenitor models for these studies. 
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APPENDIX A: PHYSICS OF 
ELECTRON-CAPTURE AND BETA-DECAY 

We are interested in the electron-capture reaction 

{Z,N) + e- + l) + v, (Al) 

and its reverse reaction, /3-decay 

{Z,N) ^ {Z +l,N -l) + e- + v, (A2) 

where Z and N are respectively the proton and neutron 
number of the nucleus. For nuclei in a dense plasma where 
the electrons are degenerate, the rates of these processes can 
depend sensitively on the density (though the electron dis¬ 
tribution function) and temperature (though the occupation 
of nuclear energy levels and the electron distribution func¬ 
tion). The neutrinos are able to free stream out of the star, 
and therefore neutrino phase-space is effectively unfilled. For 
a more thorough discussion of the physics of weak reactions 
in stellar environments, see e.g. Fuller et al. ( 1980| 19851. 

This section summarizes a simple framework for the 
rates of these weak processes. More detailed microscopic cal¬ 
culations exist in the literature such as those presented in 


Oda et al. (19941. However, those particular tables are suffi¬ 


ciently sparse that numerical considerations related to inter¬ 
polation cause us to elect to use rates calculated in the man¬ 
ner described here, rather than interpolate in tables from 
more detailed calculations. 

The rate of the electron capture or /3-decay transition 
from the i-th state of the parent nucleus to the j-th state 
of the daughter nucleus can be written as (e.g. jFuller et al.| 
19801) 


A.. = ^/(/r.T,g..), 


(A3) 


where (ft) is the comparative half-life and can be either 
measured experimentally or theoretically calculated from 
the weak-interaction nuclear matrix elements. I is a phase 
space factor which depends on the temperature T, electron 
chemical potential /r, and the energy difference Qij between 
the (i-th) parent and (j-th) daughter nuclear states. 


Qij — (/Tp — jid) + Ei — Ej 


(A4) 


where fip and jid are the chemical potentials of the nuclei. 
For a classical ideal gas, the chemical potential is 


III = m/c^ -I- kT In ( — 

. riq 


(A5) 


where mj is the rest mass, n/ is the number density, and 
Uq = {2'KmikT. Therefore, 


Qij = (Mp -Md)c^ + kT\n[Z^]EEi- E, 


rid 


(A6) 


where Mp and Md are the nuclear rest masses of the par¬ 
ent and daughter nuclei, respectively. Since \Mp — Md \ <? ~ 
5 MeV for the isotopes we consider and we restrict ourselves 
to temperatures T < 10® K (so kT < 100 keV), the term 
A:Tln^(^j is negligible in comparison and we discard it, 
leaving 


Qij — (Afp — Md) c -\- Ei — Ej 


(A7) 


Though not generally true, for the set of transitions that we 


consider, this definition means that Qij < 0 for e-capture 
and Qij > 0 for /3-decay. 

We work in the allowed approximation, which neglects 
all total lepton angular momentum {L = 0). This restricts 
us to the following transitions and corresponding selection 
rules (e.g. Commins|1973 1: Fermi transitions, where the to¬ 
tal lepton spin is S = 0 , and therefore the initial and final 
nuclear spins are equal (Ji = Jj), and Gamow-Teller transi¬ 
tions, where S = 1, and therefore Ji = Jj, Jj ± 1 (excluding 
Ji = J/ = 0). In both cases, these is no parity change: 
TTiTT/ = +1- 

At low temperature, the electron chemical potential is 
approximately the Fermi energy Ep (the first correction en¬ 
ters at order {kT/Ep)'^), and so we use the terms Fermi 
energy and chemical potential interchangeably. In the rel¬ 
ativistic limit, the chemical potential can be approximated 


where Ye = X^i ZiXijAi is the electron fraction. Zi, Xi, and 
Ai are respectively the charge, mass fraction, and atomic 
mass of the i-th species. 

The total rate of the process is the sum of the individ¬ 
ual transition rates from the i-th parent state to the j-th 
daughter state, Xij, weighted by the occupation probability 
of the i-th parent state, pi. 

Atotal ~ ^ ^ ^ ^ Aij , (-A9) 

® i 

The i-sum is over all parent states and the j-sum is over 
all daughter states. We will always assume that the nu¬ 
clear states are populated with a thermal (Boltzmann) dis¬ 
tribution. Some parent nuclei preferentially capture into ex¬ 
cited daughter states, but these excited states decay via 7 - 
ray emission with a typical timescale ~ 10“^® s. Therefore 
the level population returns to a thermal distribution on a 
timescale much shorter than the evolutionary timescales of 
interest]^ The occupation probability is 

9 I -I- 1 

W = ^f^ei^p{-PEi) , (AlO) 

where P{T), the nuclear partition function is P{T) = 
X)i(2Ji -I- and we define /3 = {kT)~^. 

The remainder of this section considers the rate of a 
single allowed {L = 0) transition in detail, and so for con¬ 
venience we drop the i,j subscripts. In the case of electron 


® There is one case in which this hierarchy of timescales is not so 
obvious. The firs t excited state o f is metastable with a half- 

life of 2 X 10“® s |Firestone|2007[ | and ground state of prefer¬ 
entially captures into this excited state. If the capture rate from 
this excited state to ^^Ne were approximately equal or greater 
than the rate of decay via 7 -ray emission, the relative state popu¬ 
lations would be effectively non-thermal. Using the parameters of 
this transition as listed in Table^ the capture timescale is approx¬ 
imately equal to the half-life at a critical density of logj^Q p Ri 10.5 
(for Ye ~ 0.5), which is safely outside of the density range that 
we consider in this work. 
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capture, the phase space factor is (e.g., Fuller et al.yi980l 


toe. - 


1 /■°° 


-Q l + exp[l3{Eo - p.)] 


G{Z,Ee)dEe , 


(All) 

along with the energy conservation relationship E^. + Q = 
Eu- The quantity G is defined as 


G{Z, Eo) = Eo) 

-C/e 


(A12) 


where E{Z,Ee) is the relativistic Coulomb barrier factor 


(Gove & Martin 19711. We make the approximation that 


the electrons are relativistic. In this limit, 

, i:, 

/ 4nEeR 


G(Z,Eo) 


exp {naZ) 


(A13) 


\ he 

where a is the fine structure constant and R is the size 


of the nucleus (Fuller et al. 19801. We are considering nu¬ 


clei with Z « 10, A ~ 20 (and so 7? « 3fm), at densi¬ 
ties such that Es ~ 5MeV. Therefore, the value of the first 

term is ^ « 0.999, with an extremely weak E^ 

and R dependence (since a^Z^ ~ 0.005). Therefore we treat 
G(Z, Ee) as a constant with a value of exp(7raZ). 

Changing to dimensionless variables e = .^^2 a-nd q = 
- 2 - 2 -, the phase space integral becomes 


he = e 


/ OO 




, o + exp[f5msC^{e - fi)] 

We rewrite this integral in simpler form as 

^7raZ 

-fee — 


de 


(A14) 


{PrrieC 


,2'i5 


[FA{n + C.)-KF3{n + 0 + eF2{r) + 0] , 

(A15) 

where we have defined the quantities v = and ( = (3q 
and 


roo 

F.M-l J 


+ exp (a; — y) 


dx 


(A16) 


is the complete Fermi integral. Evaluating the rate requires 
evaluating three complete Fermi integrals, for which efficient 
numerical routines exist (e.g. Aparicio]|1998 1 . 

In addition to the e-capture rate, we need the rate of 
energy loss via neutrinos, which for a single transition can 
be written as 


_ meC^ In 2 rr \ 
^u.ii — / « ,N J J- t Wij ) 


[ftu — — 

where J is phase space factor, defined by an integral similar 


neutrino energy: 


to equation (All I, except with an additional power of the 

G{Z,Eo)dEe . 


dec — 


1 

-eC2)6 7_„ 1 


(mec2)6 J_Q 1 -I- exp[^(F:e - u)] 

In terms of complete Fermi integrals, this is 


(A18) 


(A19) 

The total neutrino loss rate can calculated via an 
occupation-weighted average, analogous to that used to cal¬ 


culated the total rate in equation (A9l. 


In the case of /3-decay, the phase space factor is 


Ib = 


f 

«/ m, 


^2 1 -I- exp[-,d(F: - y)] 


G{Z, E)dEe , (A20) 


and energy conservation Q = Ee + E^. Following the same 
procedure as the electron capture cas^ 


T TT 

Ib = e 


\e-qy 


rde . (A21) 


h 1 + exp[—/3mec2 (e — y)] 

We can convert this integral into a sum of complete Fermi 
integrals by making use of the mathematical identity 

k 


L 




0 1 -k exp(a; - y) 


3=0 


= -Ffe(?/) - J - b) . 


(A22) 


Defining d = fimeC^, this yields the following expression: 


Ib = 


(/3meC2)5 


(PnieC'^ 


fC 

Jp= / 

J m 


[F4(C -v)- 2CF3(C -q) + CF^C - v)] - 

[F^d-n)- 

Fiid — r])x (4i9 — 2^) -k 
F 2 {d -v)x {6d^ - 6d( + 0 - 
Fi(d -v)x (4i9^ - 6d\ -k 2dC) + 
Fo{d-v) X {d‘^ - 2d\ + d\^) ] . 

(A23) 

Similarly, the factor Jp necessary to calculate the neutrino 
loss rate is 

cQ jp-d. rp'6 

- E-CE - ^G{Z,E)dEe, (A24) 

,2 f + exp[-/3iE -y)] ^ ^ > 

which can be written as 

Jp = [^5(C -v)- KFaC -v) + EFbC - v)] - 

Fi{d-v) X (5r?-3C) + 

Fzid -v)x (lOr?^ - 12dC -k 3C^) - 
F 2 {d -v)x (lOr?® - 18d^C + - E) + 

FEd -v)x - 12r3®C + - 2dC) - 

Fold -v)x (d® - 3i7\ -k 3d^C - ] ■ 

(A25) 

Given the reaction rates and neutrino energy loss rates, 
we can calculate the net heating rate of the plasma. The 
energy equation for material in the star is 

where g-ec and account for the set of weak nuclear reac¬ 
tions we are considering separately and g* includes all other 


^ It is perhaps less obvious that the assumption that the electrons 
are relativistic is justified here, given the lower integration limit. 
But because electron phase space is only empty near or above the 
Fermi energy (and p > 5 MeV at the densities of interest), the in¬ 
tegrand is only significant at the upper portion of the integration 
range where this approximation is justified. 
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heating and cooling sources such as thermal neutrino losses 
and other nuclear reactions. Under the assumption of ther¬ 
mal equilibrium, the energy released by the weak reactions 
depends only on the total reaction rate, the total neutrino 
loss rate, and the ion and electron chemical potentials. The 
energy generation rate (per capture or decay) is 

£ec = ( — + Ate)Aec “ £v,ec (A27) 

£l3 = (—+ fJ'i.z — Ate)A/3 — e,/,/3 (A28) 


where fJ^i,z, f^i,z-i are the chemical potentials of the ions 
with those charges and Hs is the chemical potential of the 
electron. Defining Qg = Qoo = {Mp — Md)c?, which implic¬ 
itly making the same assumption used to derive equation 
(Al \ , the specific energy generation rates are 


9ec — —^Eec — [(Qg + Me)Aec “ Ei/.ec] , (A29) 

<7/3 = ^£/3 = ^ [(Qg - 7te)A/J - £„,/3 ] , (A30) 


where riec and np are the number densities of the species 
undergoing capture and decay. Therefore, given a list of nu¬ 
clear levels and the (/t)-values for the transitions between 
them, we can calculate the rates of and energy generation 
rates from electron-capture and /3-decay. 

Using the above approach, it would be possible to gen¬ 
erate tables whose points are spaced sufficiently closely, 
such that interpolation would no longer incur significant er¬ 
rors. But because MESA comes with fast quadrature rou¬ 
tines to evaluate equation (A16l, it directly evaluates equa¬ 
tions ( |A15[ ), ( |A19[ ), ( |A23[ ), and ( |A25[ ) each time one of the 
weak reaction rates is needed. While this is computationally 
inefficient, the overall speed of our calculations is sufficiently 
unaffected that we chose not to optimize this. 


B1 Equation of State 


Most straightforwardly, a Coulomb term appears in the ion 
equation of state. This affects the weak reaction rates, be¬ 
cause at a fixed total pressure the electron density is lower. 
The MESA equation of state routines, which in the ther¬ 
modynamic regime of interest are based on the Helmholtz 
equation of state (Timmes & Swesty 20001, include these 
terms based on the work of Yakovlev & Shalybkov (19891. 


B2 Ion Chemical Potential 


The energy required to remove an ion of one species and 
create an ion of another species is given by the difference in 
the ion chemical potentials. Since electron-capture and beta- 
decay change the ion charge, the presence of the Coulomb 
interaction energy changes the energy difference between 
the parent and daughter nuclear states. The interaction en¬ 
ergy is negative, and so decreasing the charge of the nucleus 
(as electron-capture reactions do) requires additional energy, 
which will therefore shift the onset of electron captures to 
higher density. 

To calculate this shift, we use the excess (that is, the 
part in addition to the ideal contribution) ion chemical 
potential Pex developed in the following series of papers: 


Chabrier & Potekhin (19981; Potekhin & Chabrier (20001 


Potekhin et al. (20091. We incorporate this effect by shifting 


the value of Q, as defined in equation ( |A4[ ), by an amount 
AE = /iex.p — Mex.d- This shift. 


Q' = g -f AE 


(B2) 


then enters the calculation of the phase space factors and the 
energy generation rates. In Fig. |B1| the red dotted line la¬ 
beled “ion chemical potential”, shows the effect of including 
of these corrections. 


APPENDIX B: COULOMB CORRECTIONS 


In a dense plasma, the electrostatic interactions of the ions 
and electrons introduce corrections to the weak rates rela¬ 
tive to rates which assume a Fermi gas of electrons and an 
ideal gas of ions (as do those presented in Appendix]^. The 
leading term in the Coulomb interaction energy for ion-ion 
interactions (e.g. Shapiro fc Teukolsky[]1983| ) is 

9 

(Bl) 




10 tti 


where ai = 1 I is the inter-ionic spacing. For Z ~ 10, 
Ec « —0.2 E_f indicating that the interactions are energeti¬ 
cally importantly 

In this section, we discuss our treatment of these cor¬ 
rections, and compare our approach to previous work. Our 
treatment is most similar to that discussed in Appendix A 
of jJuodagalvis et al.| ( |2010[ ). Fig. |B1| illustrates the effects 
of including these corrections on the evolution of our fidu¬ 
cial model. The change in the evolution is similar to that 


observed by Gutierrez et al. (19961. 


The Coulomb interaction energy and the Fermi energy both 
scale oc in the relativistic limit. 


B3 Screening 

The electron density relevant to the reaction rate is not the 
average electron density, but rather the electron density at 
the position of the nucleus. [Itoh et ffi(] ( |2002l ) calculated the 
value of this screening correction using linear response the¬ 
ory. This correction can be correctly accounted for as a shift 
in the value of the electron chemical potential that enters 
the phase space factor. 

/ie = /Te -I- Vs (B3) 

However, this correction does not enter the energy gener¬ 
ation rates because it has not changed the energy cost to 
add or remove an electron from the bulk Fermi sea, which 
is the net effect of a capture or decay. In Fig. |B1[ the yellow 
dashed line labeled “electron screening”, shows the effect of 
including these corrections. 


B4 Comparison with Previous Work 


The effect discussed in §B2| which is the dominant Coulomb 
correction, has previously been included in studies of 


ONeMg cores (e.g. Gutierrez et al. 1996 Takahashi et al. 


20131. The approach taken in these studies is to include this 


effect as a shift in the electron chemical potential 
fJ-'e = fj-e — AE 


(B4) 
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Figure Bl. Illustration of the effect of Coulomb corrections on 
the evolution of the central density and temperature of the ac¬ 
creting core. (To better understand this plot the reader may first 
want to consult Fig. and the surrounding discussion.) All cal¬ 
culations include the corrections to the equation of state (§ |B1| . 
The dashed-dotted blue line shows the result with no other cor¬ 
rections. The dashed yellow line shows the effect of the inclusion 
of the screening corrections (§ |B3[ |. The dotted red line shows the 
effect of the inclusion of the corrections to the ion chemical po¬ 
tential (§ |B2[ l.The solid black line shows the result with both cor¬ 
rections included, which is the default choice for our calculations. 
The primary effect of including these corrections is a increase (of 
about 0.05 dex) in the density at which electron captures occur. 


and to use this modified electron chemical potential in the 
evaluation of the rates. This approach is conceptually incor¬ 
rect, because and Q enter the rate expression in different 
ways, as can be seen in equation (Alll. However, given a 
table of Aec(p, T), one has no ability to shift Q, so the only 
way to correct the rate is to shift the relation between /le 
and p. In the sub-threshold case (see equation]^, the most 
important term is the exponential, which is symmetric in Q 
and p, and so this approach does not lead to a substantial 
quantitative error in the rate. 

When making this correction both [Gutierrez et al. 


(19961 and Takahashi et al. (20131 follow Couch & Loumos 


(19741 and use the form of the ion free energy from Dewitt 


et al. (19731. There has been progress in calculating the free 


energy of electron-ion plasmas in the last few decades. As 
discussed in § |B2[ we use the fitting formula for the free en¬ 


ergy from Potekhin et al. (20091. In their work on electron 
capture rates in NSE material, Juodagalvis et al. ( 20101) use 


the formula quoted in Ichimarul (19931. The results ofllchl 


(19931 and Potekhin et al. (20091 agree, while the shift 


calculated following Dewitt et al. 


(19731 is approximately 30 


per cent larger in magnitude. 

The screening correction discussed in § M is not as 
widely adopted. It is included in Juodagalvis et al. (20101, 


but not in Gutierrez et al. (19961. The results of Itoh et al. 


(20021 are within approximately 10 per cent of the results 


Figure Cl. The legend shows the number of timesteps and the 
maximum number of zones used in the calculation of each of the 
runs shown in Table [C] The negligible variation between models 
indicates our results are converged. 


from the Thomas-Fermi approximation 

14 « 


(B5) 


This effect has approximately the magnitude of the differ¬ 
ence between older ion chemical potential and the one we 
adopt discussed in the preceding paragraph, but the oppo¬ 
site sign. Therefore, despite its exclusion, the net difference 


between our calculations and those of Gutierrez et al. (19961 
is small. 


APPENDIX C: CONVERGENCE 

In order to demonstrate that our results are robust, we per¬ 
form a number of tests of the spatial and temporal conver¬ 
gence of our MESA calculations. The parameters of these 
runs are shown in Table We performed runs with each 
of the spatial and temporal resolutions each separately ten 
times greater than the hducial case, as well as a run in which 
both the spatial and temporal resolutions were three times 
greater than the fiducial case. Fig. |Cl| shows that the central 
temperature evolution remains unchangecp^nd we observed 
that the variation of the result in any quantity of interest 
was negligible. 

In addition, we performed a run with a much larger 
network (203 isotopes; mesa_201 .net plus and ^'^F) and 
confirmed that our results remained unchanged. 


The small difference at the end of the “Temporal” track is an 
artifact of a difference in how the stopping condition trigger was 
tripped. 
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Run Name 

delta_lgRho_cntr_liinit 

delta_lgRho_ciitr_hard_liinit 

varcontrol_target 

mesh_delta_coeff 

Fiducial 

1 X 10-3 

3 X 10-3 

1 X 10-3 

1.0 

Temporal 

1 X 10-"* 

3 X 10-"^ 

1 X 10-* 

- 

Spatial 

- 

- 

- 

0.1 

Both 

3 X lO-'* 

1 X 10-3 

3 X 10-* 

0.3 


Table Cl. Parameters for the runs demonstrating the convergence of our results. The column names are the specific MESA controls we 
used. The controls delta_lgRho_ciitr_limit, delta_lgRho_cntr Jiard_limit, varcontrol_target control the maximal fractional change in 
physical variables, which in an implicit code like MESA controls the timestep. The control mesh_delta_coeff controls the number of 
zones used in the calculation. Fig. |Cl|shows that central temperature evolution is essentially identical for these different runs. 


APPENDIX D: TWO-ZONE WD MODELS 

This section describes the framework we use to understand 
the evolution of our MESA models after the A = 24 cap¬ 
tures have occurred. The quantitative estimates shown in 
Figs. and were made using the approach described in 
this Appendix. 

After the A = 24 electron captures have begun in the 
center of the WD, the MESA models have two zones: an in¬ 
ner “neutronized” zone in which the captures have occurred 
and Ee is lower, and an outer zone whose composition re¬ 
mains unchanged. This property of our MESA models can 
be seen in Fig. The neutronized zone is growing (in a 
Lagrangian sense) as a function of time. 

In § |D1| we write down an idealized model of a white 
dwarf with this two zone structure. In § |D2| we discuss how 
we apply these models to understand our MESA calcula¬ 
tions. 


D1 Details of the Two-Zone Model 


Following Cox (19681, we write down a simple model of a 
zero-temperature white dwarf. We assume spherical sym¬ 
metry and hydrostatic equilibrium and so solve the Poisson 
equation in spherical coordinates: 


1 d fr^dP\ ^ ^ 
dr \p dr ) ^ 


(Dl) 


We also assume the equation of state of a zero temperature, 
ideal Fermi gas, which is 


model. We also define 

_ /2A 1 

and transform to the variables 


(D6) 


r = aC, 
z = 2c'h . 

This yields the differential equation 



At the center ((" = 0), the boundary conditions are 


(D7) 

(D8) 

(D9) 


cl.(C = 0) = 1 

d-F 


dC 


At the surface (C = Cs): P - 


= 0 

C=o 

0 and so 2 


(DIO) 

(Dll) 


1, meaning 


1>(C = 0) = - . 


(D12) 


Now, we divide the white dwarf into two zones which 
have different values of W. By assuming a piecewise con¬ 
stant form for Te, we can continue to solve equation (D9l 
throughout the whole white dwarf. Specihcally, we use 


W = 


Ee,0 

Ee.n 


if 2 < 2n 

if 2 > 2„ 


(D13) 


The transition between the two zones occurs at C,n such that 
‘l’(Cn) = Znjza- The following physical conditions must be 
satished at this interface 


P = Af{x) 
B 


1/3 


where 


f{x) = x\/x^ + 1 (2x^ — 3) -I- 3sinh '^(x) 


(D2) 

(D3) 

(D4) 


A = 


3h3 


6.0 X 10^^ dynes cm ^ 


Combining this equation of state with equation (|D1|) gives 

. 1/2 


J__d 

r2 fij. 


2 d , 2 , 

+ 1 ) 


TtGB^ 3 

■ 2Ay2 ® 


(D5) 


where we have made the assumption that dY^/dr — 0. 

In order to non-dimensionalize these equations, define 
2^ = -I- 1 and let Zc be the value of 2 at center of the 


P+=P- 


1(W 
p dr 


f GMr 


(D14) 

(D15) 


Note that the continuity of P implies the continuity of x, and 
hence 2, even though Ye is discontinuous. The dimensionless 
equivalents of these conditions are 


$(c = a“) = ‘&(c = Cn) 




d<I>\ 


= Ye 


d-FX 


(D16) 

(D17) 


C=C4 


Constructing a two-zone model is now simple. Specify 
the three parameters for the equation of state: Yep, Ye,n, Zn- 
Select a central density (which sets the value of 2c) and then 
integrate the ODE observing the boundary and jump condi¬ 
tions. The solution gives the structure of a single two-zone 
model. A one-parameter family of models can be constructed 
by varying 2c, which in turn varies the properties (e.g., mass, 
radius) of the model. 
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D2 Applications of the Two-Zone Model 


Our MESA models are in hydrostatic equilibrium: their 
evolution is occurring on timescales much longer than the 
dynamical time. We use the two-zone model to find ap¬ 
proximate sequences of hydrostatic models along which the 
MESA models evolve. This gives us insight into the processes 
that control the timescale of the evolution. 


The piecewise equation of state given in equation ( |D13[ ) 
can be used to represent the A = 24 captures: setting 
log^gP„ = 9.6, with Ve.o = 0.5 and Ve.n = - XMg/12 

corresponds to instantaneous neutronization of all available 
^■^Mg at densities above the threshold density. 


Fig. |D1| shows the schematic evolution of models with 
Aivig = 0.05. The black line is the family of two-zone hydro¬ 
static models. This family of models is generated by vary¬ 
ing Pc (the central pressure). Because of the discontinuity 
in Yc, the continuous variation in Pc gives a discontinuous 
variation in pc- (In Fig. |D1| the density jump is hidden by 


point 2, but the jumps are apparent in Fig. D2 ) The grey 
line shows the family of models without neutronization. For 
central densities less than p„ (e.g., point 1) the two families 
are equivalent. At point 2, the central density crosses the 
threshold density and the families diverge. At point 3, the 
model has a substantial low-Ye core and has a much higher 
central density at fixed mass relative to the models without 
neutronization. 


The numbered points in Fig. |D1| represent a temporal 
sequence of models. The time evolution is driven by the in¬ 
crease in M set by accretion. For a given value of M, the 
timescale for changes in pc is given by equation where the 
value of d\npc/d\nM comes from the sequence of hydro¬ 
static models. In Fig. the black dashed line is calculated 
using the model without neutronization, while the black dot¬ 
ted line is calculated using the two zone model. This latter 
line does an excellent job of quantitatively describing the 
more rapid contraction of the MESA model following the 
A = 24 captures. 

Fig. |D2| shows the schematic evolution of models with 
Xivig = 0.15. The black line is the family of two-zone hydro¬ 
static models. The grey line shows the family of models with¬ 
out neutronization. For central densities less than (e.g., 
point 1) the two families are equivalent. At point 2, neu¬ 
tronization begins to have an effect. At point 3, the model 
passes the maximum mass possible for the family of models 
with logjQ pn = 9.60. Up to this point, as in the case with 
Aivig = 0.05, the points 1-3 represent a temporal sequence 
of models whose time evolution is driven by increasing M. 


With the equation of state held fixed, point 3 would 
mark the onset of dynamical instability. However, the char¬ 
acteristic timescale of the electron captures is longer than 
the dynamical time. Therefore it is not physically possible 
for pn to remain fixed. Only material at densities where the 
electron capture timescale is shorter than the evolutionary 
timescale is able to completely neutronize. 

For the large Xyig, models like that in Fig. |D2[ the evo¬ 
lutionary timescale becomes sufficiently short that there is 
no longer time for a significant amount of mass to accrete. 
Therefore, the evolution switches to a sequence of models 
with constant M, as indicated by the black dashed line and 
points 4-6 in Fig. |D2[ The fixed value of M is the maximum 



logio Pc [g cm 3] 


Figure Dl. The black line shows the sequence of two-zone hydro¬ 
static models with electron captures above logjQ p„ = 9.6 with 
the change in Yc corresponding to XMg = 0.05. The grey line 
shows a zero-temperature model without neutronization. With 
neutronization taken into account, the central density increases 
more rapidly with increasing mass. The numbered points indicate 
a temporal series of models; specific points are discussed in more 
detail in the text. 


mass for a model with the initial value of logj^Q pn = 9.60 
(i.e., point 3). 

The evolution along the temporal sequence of points 3- 
6 is limited by the electron capture rates. To describe this 
quantitatively, for a given pc, we find the value of pn that 
corresponds to the hydrostatic model with a given M. Then 
we calculate the neutronization timescale (equation |23| l cor¬ 
responding to this value of pn- This gives the black dash 
dotted line shown in Fig. which does a good job of re¬ 
producing the evolution observed in the MESA model. For 
simplicity, we used the electron capture rates at fixed tem¬ 
perature (logjQT = 8.6) to calculate Aec. 

The qualitatively different evolution experienced by the 
Aivig = 0.05 and Amg = 0.15 models is due to the presence 
of a maximum mass in the families of two-zone models at 
a central density less than the critical density for ^°Ne cap¬ 
tures. Calculating the central density at which this maxi¬ 
mum occurs for each value of Ajug gives the critical curve 
shown as a black dashed line in Fig. 
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logio Pc [g cm 3] 


Figure D2. The black line shows the sequence of two-zone hydro¬ 
static models corresponding to X]vig = 0.15. The grey line shows 
a zero-temperature model without neutronization. The colored 
lines show models with neutronization at different densities. The 
gaps between the colored lines and the grey line are the discon¬ 
tinuities in pc caused by the discontinuity in Ye\ the sequence in 
continuous in Pc- The numbered points indicate a temporal series 
of models; specific points are discussed in more detail in the text. 
The black line changes from solid to dashed when the sequence 
of models changes from being defined by constant pn to constant 
M. 
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